function residuals = xK_iterative_residual( xKs,  kappap,  par  ) 
% Usage: xK_iterative_residual( xKs,  kappap,  par  ) 
% 
% Note:  Input xKs and kappap are either vectors of the same length, or kappap could be a scalar, 
% in which case the algorithm construct an array of kappap with the same length of xK. 

lambda=par.lambda;  alpha=par.alpha;  beta=par.beta;   varepsilon=par.varepsilon;
sigmaK=par.sigmaK; AH=par.AH; AL=par.AL; 
yK_fun=par.yK_fun;  sigmap_fun=par.sigmap_fun;  
p = par.p_ind ;    w=par.w_ind;
xK_max = par.xK_max;    xK_min = par.xK_min ;

if(numel(kappap)~=numel(xKs) )
    kappap = kappap(1)*ones( numel(xKs), 1 );
end

yKs = yK_fun(xKs);
sigmap = sigmap_fun(xKs);

% Previous version
% residuals = ( AH-AL )/p - ( sigmaK + sigmap  ).^2.*( xKs - yK ) - lambda.*( 1./( (1+alpha*beta)./(kappap+alpha*beta)  - xKs )  -  kappap./( 1- yK.*kappap )  ) ;
delta_x = beta*max(xKs-1,0);
indicator = xKs -1 >0 ; 
kappa_fs = alpha*delta_x*w/(1-w);
kappad = varepsilon;

residuals = ( AH-AL )/p - ( sigmaK + sigmap  ).^2.*( xKs - yKs ) ...
               -  lambda*( (kappap+ indicator.*alpha*beta)./( 1-xKs.*kappap-alpha.*delta_x ) ) ...
              +  lambda*( ( kappap - kappad ) ./ (  1-yKs.*kappap-(1-yKs).*kappad + kappa_fs )  )  ...
              +  10^10*max(xKs-xK_max,0) + 10^10*max(xK_min - xKs,0);


